clear;clc;
tic;

%% 载入图案数据
% 共有十组数据
path = pwd; % 获取当前m文件夹路径
filepath = [path,'\generate_shape\generate_top_bottom\'];
groupfilepath1 = ['top_arc_bottom_arc'];
groupfilepath2 = ['top_arc_bottom_rectangle'];
groupfilepath3 = ['top_Hshape_bottom_arc'];
groupfilepath4 = ['top_Hshape_bottom_rectangle'];
groupfilepath5 = ['top_rectangle_bottom_arc'];
groupfilepath6 = ['top_rectangle_bottom_Hshape'];
groupfilepath7 = ['top_rectangle_bottom_rectangle'];
groupfilepath8 = ['top_rectangle_bottom_Xshape'];
groupfilepath9 = ['top_Xshape_bottom_arc'];
groupfilepath10 = ['top_Xshape_bottom_rectangle'];
finalpath = groupfilepath2;
topfile = [filepath,finalpath,'\top.mat'];
bottomfile = [filepath,finalpath,'\bottom.mat'];

% top = importdata(topfile); % 下层图案
% bottom = importdata(bottomfile); % 上层图案
top = importdata('top.mat');
bottom = importdata('bottom.mat');
for n=1:8
    %% 调用CSTStudio并初始化
    a = top(:,:,n);
    A = bottom(:,:,n);
    cst = actxserver('CSTStudio.application');
    mws = invoke(cst,'NewMWS');
    app = invoke(mws,'GetApplicationName');
    ver = invoke(mws,'GetApplicationVersion');
    invoke(mws,'FileNew');
    % filename = ['\metamaterials\metamaterial',num2str(n),'.cst']; % 新建的CST文件名字
    % fullname = [path filename];
    % invoke(mws,'SaveAs',fullname,'True'); % 保存到目前为止的结果
    invoke(mws,'DeleteResults'); %删除之前的结果，防止程序停止运行

    %% 仿真参数设置
    Px = 64; %结构X向边长
    Py = 64; %结构Y向边长
    h = 2; %基地层金属厚度
    h1 = 0.5; % 顶层金属层厚度
    h2 = 0.5; %底层金属层厚度
    l1 = 35; %顶层金属边长
    l2 = 50; %底层金属边长
    w = 5; %金属宽度
    t1 = 5; %顶层介质层厚度
    t2 = 5; %底层介质层厚度
    Frq = [2.5,3.5]; % 工作频率
    little_length = 1; % 图案小块边长

    %% 单元格式设置
    units = invoke(mws,'Units');
    invoke(units,'Geometry','um');
    invoke(units,'Frequency','THz');
    invoke(units,'Time','ns')
    invoke(units,'TemperatureUnit','kelvin');
    release(units);

    %% 频率范围设置
    solver = invoke(mws,'Solver');
    invoke(solver,'FrequencyRange',Frq(1),Frq(2));
    release(solver);

    %% 材料设置
    % 1. 背景材料设置
    background = invoke(mws,'Background');
    invoke(background,'ResetBackground');
    invoke(background,'Type','Normal');
    release(background);

    % 2. 材料1Gold设置
    material = invoke(mws,'Material');
    material_gold = 'Gold';
    sigma = 4.561e+007;
    invoke(material,'Reset');
    invoke(material,'Name',material_gold);
    invoke(material,'FrqType','all');
    invoke(material,'Type','Lossy metal');
    invoke(material,'Sigma',sigma);
    invoke(material,'Mu','1.0');
    invoke(material,'Rho','19320');
    invoke(material,'ThermalConductivity','314');
    invoke(material,'HeatCapacity','0.13');
    invoke(material,'YoungsModulus','78');
    invoke(material,'PoissonsRatio','0.42');
    invoke(material,'ThermalExpansionRate','14');
    invoke(material,'Colour','1','1','0') % 设置颜色为黄色
    invoke(material,'Create');
    invoke(material,'Delete','Vacuum'); % 删除原本有的两个材料
    invoke(material,'Delete','PEC')
    release(material);

    % 3. 材料Dielectric设置
    material = invoke(mws,'Material');
    material_dielectric = 'Dielectric_material';
    Mu = 2.0;
    invoke(material,'Reset');
    invoke(material,'Name',material_dielectric);
    invoke(material,'FrqType','all');
    invoke(material,'Type','Normal');
    invoke(material,'Epsilon','1.0');
    invoke(material,'Mu','2.0');
    invoke(material,'Create');
    release(material);

    % 4. 建立新的部分
    component = invoke(mws,'Component');
    invoke(component,'New','Substrate');
    invoke(component,'New','Plasmonic');
    invoke(component,'New','Dielectric');

    %% 建立模型
    %% 1. 建立底层金属模型
    Name = 'Substrate';
    component1 = 'Substrate';
    material = 'Gold';
    brick = invoke(mws,'Brick');
    invoke(brick,'Reset');
    invoke(brick,'Name',Name);
    invoke(brick,'Component',component1);
    invoke(brick,'Material',material);
    invoke(brick,'Xrange',0,Px);
    invoke(brick,'Yrange',0,Py);
    invoke(brick,'Zrange',-h,0);
    invoke(brick,'Create');
    release(brick);

    %% 2. 建立下层介质层模型
    Name = 'bottom Dielectric';
    component2 = 'Dielectric';
    material = 'Dielectric_material';
    brick = invoke(mws,'Brick');
    invoke(brick,'Reset');
    invoke(brick,'Name',Name);
    invoke(brick,'Component',component2);
    invoke(brick,'Material',material);
    invoke(brick,'Xrange',0,Px);
    invoke(brick,'Yrange',0,Py);
    invoke(brick,'Zrange',0,t2);
    invoke(brick,'Create');
    release(brick);

    %% 3.1. 建立下层金属层
    k = 1; % 用于统计小块的个数
    for i = 1:64
        for j = 1:64
            if A(i,j) ==1
                Name = ['Ring', num2str(k)];
                Name_paste = ['Ring',num2str(k),'_1'];
                brick = invoke(mws,'Brick');
                invoke(brick,'Reset');
                component3 = 'Plasmonic';
                material = 'Gold';
                invoke(brick,'Name',Name);
                invoke(brick,'Component',component3);
                invoke(brick,'Material',material);
                Xrange = [little_length*(i-1) little_length*(i)];
                Yrange = [little_length*(j-1) little_length*(j)];
                Zrange = [t2 t2+h2];
                invoke(brick,'Xrange',Xrange(1),Xrange(2));
                invoke(brick,'Yrange',Yrange(1),Yrange(2));
                invoke(brick,'Zrange',Zrange(1),Zrange(2));
                invoke(brick,'Create');
                release(brick);
                k = k+1;
            end
        end
    end

    % 3.2. 合并下层金属
    for i=1:k-2
        solid = invoke(mws,'Solid');
        invoke(solid,'Add',['Plasmonic:Ring','1'],['Plasmonic:Ring', num2str(i+1)]);
        release(solid);
    end

    %% 4.1. 创建上层介质层
    Name = 'top Dielectric';
    component2 = 'Dielectric';
    material = 'Dielectric_material';
    brick = invoke(mws,'Brick');
    invoke(brick,'Reset');
    invoke(brick,'Name',Name);
    invoke(brick,'Component',component2);
    invoke(brick,'Material',material);
    invoke(brick,'Xrange',0,Px);
    invoke(brick,'Yrange',0,Py);
    invoke(brick,'Zrange',t2,t1+t2);
    invoke(brick,'Create');
    release(brick);
    
    % 4.2 在上层介质层中创建相同的下层金属层并除去
    Transform = invoke(mws,'Transform');
    invoke(Transform,'Reset');
    invoke(Transform,'Name','Plasmonic:Ring1');
    invoke(Transform,'Vector','0','0','0');
    invoke(Transform,'UsePickedPoints','False');
    invoke(Transform,'InvertPickedPoints','False');
    invoke(Transform,'MultipleObjects','True');
    invoke(Transform,'GroupObjects','False');
    invoke(Transform,'Repetitions','1');
    invoke(Transform,'MultipleSelection','False')
    invoke(Transform,'Destination','');
    invoke(Transform,'Material','');
    invoke(Transform,'Transform','Shape','Translate');
    release(Transform);
    
    Solid = invoke(mws,'Solid');
    invoke(Solid,'ChangeComponent','Plasmonic:Ring1','Dielectric');
    invoke(Solid,'Subtract','Dielectric:top Dielectric','Dielectric:Ring1');
    release(Solid);

    %% 5.1. 创建顶层金属层
    clear k;
    k = 1;
    for i = 1:64
        for j = 1:64
            if a(i,j) ==1
                brick = invoke(mws,'Brick');
                invoke(brick,'Reset');
                Name = ['top metallic', num2str(k)];
                component3 = 'Plasmonic';
                material = 'Gold';
                invoke(brick,'Name',Name);
                invoke(brick,'Component',component3);
                invoke(brick,'Material',material);
                k = k+1;
                Xrange = [little_length*(i-1) little_length*(i)];
                Yrange = [little_length*(j-1) little_length*(j)];
                Zrange = [t1+t2 t1+t2+h1];
                invoke(brick,'Xrange',Xrange(1),Xrange(2));
                invoke(brick,'Yrange',Yrange(1),Yrange(2));
                invoke(brick,'Zrange',Zrange(1),Zrange(2));
                invoke(brick,'Create');
                release(brick);
            end
        end
    end

    % 5.2. 合并顶层金属
    for i=1:k-2
        solid = invoke(mws,'Solid');
        invoke(solid,'Add',['Plasmonic:top metallic','1'],['Plasmonic:top metallic', num2str(i+1)]);
        release(solid);
    end
    clear k

    %% 定义边界条件及端口
    FloquetPort = invoke(mws,'FloquetPort');
    invoke(FloquetPort,'Reset');
    invoke(FloquetPort,'Port','Zmin');
    invoke(FloquetPort,'SetNumberOfModesConsidered','2');
    invoke(FloquetPort,'SetUseCircularPolarization','True');
    invoke(FloquetPort,'Port','Zmax');
    invoke(FloquetPort,'SetNumberOfModesConsidered','2');
    invoke(FloquetPort,'SetUseCircularPolarization','True');
    release(FloquetPort);

    Boundary = invoke(mws,'Boundary');
    invoke(Boundary,'Xmin','unit cell');
    invoke(Boundary,'Xmax','unit cell');
    invoke(Boundary,'Ymin','unit cell');
    invoke(Boundary,'Ymax','unit cell');
    invoke(Boundary,'Zmin','expanded open');
    invoke(Boundary,'Zmax','expanded open');
    release(Boundary);

    %% 设置求解器求解
    invoke(mws,'ChangeSolverType','HF Frequency Domain')
    FDSolver = invoke(mws,'FDSolver');
    invoke(FDSolver,'Reset');
    invoke(FDSolver,'AccuracyTet','1e-3');
    invoke(FDSolver,'AccuracySrf','1e-3');
    invoke(FDSolver,'LimitIterations','True');
    invoke(FDSolver,'MaxIterations','30');
    invoke(FDSolver,'Type','Iterative');
    invoke(FDSolver,'SetNumberOfResultDataSamples','1001');
    invoke(FDSolver,'SetResultDataSamplingMode','Frequency (linear)');
    invoke(FDSolver,'Start');

    %% 后处理
    CST_Post = CST_MicrowaveStudio();
    [freq,sparam,stype] = CST_Post.getSParameters;
    % 得到的S参数依次为：
    % Szmax1,zmax1    Szmax2,zmax1    Szmin1,zmax1    Szmin2,max1
    % Szmax1,zmax2    Szmax2,zmax2    Szmin1,zmax2    Szmin2,zmax2
    % Szmax1,zmin1     Szmax2,zmin1     Szmin1,zmin1    Szmin2,zmin1
    % Szmax1,zmin2     Szmax2,zmin2    Szmin1,zmin2    Szmin2,zmin2
    % Szmax1,zmax1 == Szmax2,zmax2
    % 取CD = Szmax2,zmax1-Szmax1,zmax2
    % frequency = freq(:,1); % 获取频率点
    % save('frequency.mat','frequency');
    real = real(sparam); % 获取实部
    imag = imag(sparam); % 获取虚部
    linear = sqrt(real.*real+imag.*imag); % 计算linear值
    Szmax1zmax1 = linear(:,1);
    Szmax2zmax1 = linear(:,2);
    Szmax1zmax2 = linear(:,5);
    CD = Szmax2zmax1.*Szmax2zmax1-Szmax1zmax2.*Szmax1zmax2;
    real_name = [path, '\result\real\real_',num2str(n),'.mat'];
    imag_name = [path, '\result\imag\imag_',num2str(n),'.mat'];
    Szmax1zmax1_name = [path, '\result\S-Parameters\Szmax1zmax1_',num2str(n),'.mat'];
    Szmax2zmax1_name = [path, '\result\S-Parameters\Szmax2zmax1_',num2str(n),'.mat'];
    Szmax1zmax2_name = [path, '\result\S-Parameters\Szmax1zmax2_',num2str(n),'.mat'];
    CD_name = [path, '\result\CD\CD_',num2str(n),'.mat'];
    save(real_name,'real');
    save(imag_name,'imag');
    save(Szmax1zmax1_name,'Szmax1zmax1');
    save(Szmax2zmax1_name,'Szmax2zmax1');
    save(Szmax1zmax2_name,'Szmax1zmax2');
    save(CD_name,'CD')
    clear freq sparam stype real imag linear Szmax1zmax1 Szmax2zmax1 Szmax1zmax2 CD 
    
    %% 保存并退出
    filename = ['\metamaterials\',finalpath,'\metamaterial',num2str(n),'.cst']; % 新建的CST文件名字
    fullname = [path filename];
    invoke(mws,'SaveAs',fullname,'True'); % 保存到目前为止的结果
    invoke(mws,'Quit');
end
toc